Stochastic Schroedinger equation from optimal observable evolution 
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In this article, we consider a set of trial wave- functions denoted by \Q) and an associated set 
of operators Ac, which generate transformations connecting those trial states. Using variational 
principles, we show that we can always obtain a quantum Monte-Carlo method where the quantum 
evolution of a system is replaced by jumps between density matrices of the form D — \Qa) {Qb\, and 
where the average evolutions of the moments of Ac up to a given order k, i.e. (Aai), {AaiAa2),---, 
{Aai ■ ■ ■ Ao,^}, are constrained to follow the exact Ehrenfest evolution at each time step along each 
stochastic trajectory. Then, a set of more and more elaborated stochastic approximations of a 
quantum problem is obtained which approach the exact solution when more and more constraints 
are imposed, i.e. when k increases. The Monte-Carlo process might even become exact if the 
Hamiltonian H applied on the trial state can be written as a polynomial of Ac- The formalism 
makes a natural connection between quantum jumps in Hilbert space and phase-space dynamics. 
We show that the derivation of stochastic Schroedinger equations can be greatly simplified by taking 
advantage of the existence of this hierarchy of approximations and its connection to the Ehrenfest 
theorem. Several examples are illustrated: the free wave-packet expansion, the Kerr oscillator, a 
generalized version of the Kerr oscillator, as well as interacting bosons or fermions. 
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I. INTRODUCTION 

Starting from the pioneering works of Feynman on 
path integrals [ij, large efforts have been devoted to the 
possibility to replace a quantum problem by a stochastic 
process IE Hi IE IE Q ■ One of the greatest interest of 
stochastic formulations can be summarized as follow: let 
us consider a quantum system where the number of de- 
grees of freedom to follow is too high to authorize an ex- 
act treatment. In some cases, it is possible to recover the 
exact solution by averaging over an ensemble of trajecto- 
ries where the number of degrees of freedom to consider 
along each path is significantly reduced. Then the exact 
problem can be recast into a set of problems that can be 
carried out. Such a stochastic formulation is for instance 
at the heart of quantum Montc-Carlo techniques in par- 
ticular when considering interacting systems |E IE • 

Essentially two strategies exist to introduce stochastic 
formulations. The first technique consists of selecting a 
limited number of degrees of freedom associated to spe- 
cific observable, which contains already a large amount of 
information on the system. Then, the system evolution is 
recovered by averaging over stochastic evolution of these 
observable. This technique will be called here as 'Phase- 
Space technique' |E|- New developments along this line 
have been proposed to treat interacting fermionic systems 
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Alternatively, several works have shown that one can 
take advantage of the Stochastic Schroe ding er Equation 

fE) a£proach_p, IlLI, MlIilllllllllllEliE 
120I I2II to simulate exactly quantum systems 
123. 124L I25L l2a | . In that case, a specific class of trial 
wave-functions is selected and the exact Schroedinger 
equation is mapped into a stochastic process between 
trial states. Again, trial wave- functions are chosen to 



be much simpler than the exact wave-function, although 
complex enough to incorporate already some important 
features of the problem. SSE have generally an equiv- 
alent stochastic evolution in phase-space: for instance, 
the stochastic process proposed in ref. p3 | can be equiv- 
alently replaced by a stochastic process for the one-body 
density evolution |23| associated with N-body densities 
D = |<I>i) ($2! where both states corresponds to trial 
wave- functions, i.e. Hartree-Fock states. 

In this work, we would like to address more system- 
atically the possibility to obtain SSE for a given quan- 
tum problem using the phase-space evolution. Following 
ref. [23, we consider trial state vectors denoted by \Q), 
where Q = {qa}a=i n (with TV eventually infinite) corre- 
spond to a set of variables which completely determines 
the state. The ensemble of states obtained by taking 
all possible values of Q is denoted by Sq. The second 
important hypothesis is that we assume the existence of 
an ensemble of operators {Aa}^^i j^, that generate local 
transformations between the state of Sq, i.e. one state 
\Q) transforms into another state \Q + SQ) ofSg through 
the transformation 



|Q + <5Q) =eS='«°^° IQ), 



(1) 



At the deterministic level, variational principles |29l |30j 
provide a systematic way to replace the exact problem by 
an approximate dynamics in restricted space Sq . A first 
example of extension of variational methods to obtain 
stochastic Schroedinger equation has been given in ref. 

In this article we propose an alternative procedure to 
obtain SSE from variational principles. Some aspects re- 
lated to standard variational methods are first recalled. 
Schroedinger equations obtained in this way are inti- 
mately connected to the phase-space evolution and can 
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be regarded as the best projected approximation for the 
dynamics within the subspace of observable (Aa). We 
show that variational techniques can be helpful to ob- 
tain stochastic formulations of a given problem which 
has a natural interpretation in terms of observable evo- 
lution. Starting from variational methods, we prove the 
existence of a hierarchy of approximations of a quantum 
problem using SSE methods where the expectation val- 
ues of the moments of the up to a given order fc, i.e. 
(^ai), {Aa.Aa^),..., {A^^ ' ' - ^aj, foUow in average the 
exact Ehrenfest evolution at each time step along each 
stochastic trajectory. Therefore, this generalization of 
the deterministic case can be regarded as the best ap- 
proximation of the dynamics on an enlarged subspace of 
observable. 

In the second part of this article, the existence of such 
a hierarchy and its connection to the Ehrenfest evolution 
are used to deduce SSE in several illustrative examples. 



II. GENERAL REMARKS ON VARIATIONAL 
PRINCIPLES 

Variational principles are powerful tools to provide ap- 
proximate solutions for the static or dynamical properties 
of a system when few deg rees of freedom are supposed to 
dominate m US IH Hi HI Hi ■ The selection of few 
relevant degrees of freedoms generally strongly reduces 
the complexity of the considered problem. Before intro- 
ducing stochastic mechanics, we consider standard varia- 
tional principles and recall specific aspects which will be 
helpful in the following discussion. We assume that the 
system is given at time to by |5'(io)) = \Q) where |Q) 
belongs to Sq. The wave-function, denoted by 
which evolves according to a given Hamiltonian H has 
a priori no reason to remain in the subspace Sq. How- 
ever, a good approximation of the dynamics restricted 
to wave-functions of Sq is obtained by minimizing the 
actionjil 113 

5= r ds{Q\indt~H\Q), (2) 

J to 

with respect to the variation of Q and with the additional 
boundary condition |(5<i>(to)) = and = 0. The 

variation with respect to {5Q\ leads then to the equation 

{5Q\ind'i~H\Q)^Q, (3) 

while the variation with respect to \5Q) gives the hermi- 
tian conjugate equation 

{Q\ind^ + H\6Q)^G. (4) 

Here 9^ and means that the derivative is acting re- 
spectively to the left and right hand side. 



A. observable evolution 

Due to relation J^l, the above Schroedinger equa- 
tion corresponds to specific evolutions of the expecta- 
tion values {Aa) of the operators A^. Using \5Q) = 
^^5qaAa\Q), variations with respect to \8Q) are re- 
placed by a set of independent variations {5qa}fy.^i ^, 
leading to: 

ih{dQ\A^\Q)^^dt{Q\HA^\Q), (5) 
while variation with respect to 5q'^ gives 

ih {Q\A^\ dQ) ^dt(Q\A^H\Q) . (6) 
Combining the two equations gives 

ih^^ ^ {[A^,H]), (7) 

which is nothing but the exact Ehrenfest equation of mo- 
tion for the set of operators A^. Therefore, starting from 
a density D{to) = \Q) {Q\, for one time step the evo- 
lution of < > corresponds to the exact evolution 
although the state is constrained to remain in the Sq 
space. This property is however a priori only true for 
observable which arc linear combination of Aa. 



B. Links with effective Hamiltonian and projected 
dynamics 

The evolution deduced from the variational method 
can equivalently be interpreted as an effective Hamilto- 
nian dynamics. Indeed, we have 

\dQ) = ^ dqaAa \Q) = ^i/i \Q) , (8) 

a 

where we have introduced the effective Hamiltonian Hi. 
This Hamiltonian differs from H unless H applied on the 
trial wave- function is itself a linear combination of Aa. 
The expression of Hi can be obtained using the evolution 
of dqa . From eq. ((BJ we obtain 

in J2 dqa. {ApAa) = dt {ApH) . (9) 

a 

Writing {Q\AfjAa\Q) = Mpa, dqa can formally be ob- 
tained by inverting the last expression into 

'^9" = iE^a/3'(^? 1^/3^1 (10) 

Plugging last equation into JHJl leads to Hi = ViH, where 
Vi is defined as 

ri=J2Ac.\Q)Mj{Q\A0, (11) 

a/3 
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and plays the role of a "projector" onto the subspace of 
operators A^. In particular, for any observable B, the 
operator B' = {1 - 'Pi)B fulfills {A^B') = 0. Similarly, 
if we denote by B" = B{1 ~ Vi), we have {B" A^) = 0. 
Writing the full Hamiltonian as 

H = ViH+{l-Vi)H, 

we see that the effective evolution of \dQ) given by eq. 
where Hi is obtained from the minimization of the 
action can be interpreted as the approximate evolution 
when the part of the Hamiltonian that do not contribute 
to the {Aa) is neglected. As a matter of fact, the second 
part is responsible from the deviation of the projected 
dynamics from the exact one. In summary, we have re- 
called here that the variational dynamics can be regarded 
as a projected dynamics into the subspace of relevant ob- 
servable [3^ I4Q] . In the following, we show that aspects 
of variational principles presented here can be adapted 
to obtain Monte-Carlo methods. 



III. STOCHASTIC QUANTUM MECHANICS 
FROM VARIATIONAL PRINCIPLES 

Several recent studies have shown that the exact evolu- 
tion of complex quantum systems can be replaced by a set 
of stochastic evolutions of simpler trial wave- functions. 
In this context, at a given intermediate time, the exact 
densi ty i s rep laced by the average over " densities" of the 

form mm El mil 

D{t)^\Qa){Qb\, (12) 

where both states belongs to Sq. In all cases, the main 
advantage of these reformulations was to prove that the 
stochastic methods can be performed imposing that both 
states remain in Sq. As in previous section, this condi- 
tion implies that we restrict variations of each state to 

\Qa + 5Qa) - e^='^°'^°|Qa), (13) 

\Qb + SQb) = e^=''°'^°|g6), (14) 

where now Sq^a^ and Sqa^ may also contain a stochastic 
part. 

The aim of the present section is to show that, given 
a class of trial states, a hierarchy of Monte-Carlo for- 
mulations of a quantum problem can be systematically 
obtained using variational methods. The description of 
the system can be gradually improved by introducing a 
set of noises, written as 

sql^^ = hi + '^d" + ^d" + • • • I,., 

. [b]* r b* , f [2] , c [3] , y^^l 

where the second, third... terms represent stochastic 
variables added on top of the deterministic contribution. 
Those are optimized to not only insure that the average 



evolution of (Aa) matches the exact evolution at each 
time step but also that the average evolutions of higher 
moments (AaAp), (AaApA^),... follow the exact Ehren- 
fcst dynamics. 

A. Step 1: deterministic evolution 

Assuming first that stochastic contributions s}a and 
77q' are neglected in eq. (|15|l . we show how variational 
principles described previously can be used for mixed 
densities given by eq. (|12|l . It is worth noticing that 
variational principle have also been proposed to estimate 
transition amplitudes (see also discussion in ,41j). In 
that case, different states are used in the left and right 
hand side of the action. This situation is similar to the 
case we are considering. We are interested here in the 
short time evolution of the system, therefore we disregard 
the time integral in equation ^ and consider directly the 
action 

S = Tr {{ihd'^ ~ ihd^ -H}D), (16) 

where Tr(.) stands for the usual trace. Starting from 
the above action, one can generalize the different aspects 
discussed in sectionmto the case of densities formed with 
couples of trial states. For instance, the minimization 
with respect to the variations {5Qb\ and \6Qa) leads to 
the two conditions 

ih{Qb\Aa\dQa) = {Qb\AaH\Qa), 

(17) 

in{dQb\Aa\Qa) - {Qb\HAa\Qa), 

from which we deduce that 

ihj^{Aa)^{[Aa,H]), (18) 

where {A^) = {Qb\Aa\Qa)- Therefore, the minimiza- 
tion of the action again insures that the exact Ehrenfest 
evolution is followed by the A^ observable over one time 
step. Similarly, the evolution of both \Qa) and \Qb) are 
given by 

\dQa) = Ec.dq?.Aa\Qa) =§VlH\Qa) 

{dQb\ = {Qb\EadqVAa^-§{Qb\HVi 
where Vi now reads 

Pi = ^ A„ \Qa) M^^ {Qb\ Ap. (19) 

In opposite to previous section, Vi cannot be interpreted 
as a projector onto the space of observable due to the 
fact that Map = {Qb \ AaAp \ Qa) is not anymore a met- 
ric for that space. However, we still have the property 
that the total Hamiltonian can be split into two parts, 
one which corresponds to the effective dynamic solution 
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of the minimization, and one which is neglected, which 
identifies the irrelevant part for the (Aq.) evolutions. It 
is worth noticing that the variational approximation pre- 
sented here for mixed densities has been already used 
without justification in ref. (26iJ to optimize stochastic 
trajectories. 



with respect to the variations of different parameters, i.e. 
(5<;°, i5q^ , ScUap and Saap- In the following, a formal solu- 
tion of the minimization procedure is obtained. We show 
that variational principles applied to stochastic process 
generalize the deterministic case by imposing that not 
only that expectation values (Aa) but also the second 
moments (AaAp), follow the exact Ehrenfest evolution. 



B. Step 2 : Introduction of gaussian stochastic 
processes 

In this section, we show that the description of the dy- 
namics can be further improved by introducing diffusion 
in the Hilbert space Sq . We consider that the evolutions 



of ' and qa^ now read 



de\ 



where c?^q ' and dr]a^ correspond to two sets of stochastic 
gaussian variables (following Ito rules of stochastic cal- 
culus with mean values equal to zero and variances 
verifying 



diMf = 
d^^%f = 



dcTap, 





(20) 
(21) 
(22) 



We assume that duap and daap are proportional to dt. 
Note that equation (|22|l reflects that the stochastic con- 
tributions to the evolutions of \Qa) and {Qb\ are inde- 
pendent. The advantage of introducing the Monte-Carlo 
method can be seen in the average evolutions of the 
states. Keeping only linear terms in dt in eq. (|14|l gives 
for instance 



\dQa) 



dq^A, 



- ^ duJa/J {Ac 
a</3 



Ap + ApA^) ) m) 



In the previous section, we have recalled that using trial 
state leads to an approximate treatment of the dynam- 
ics associated to effective Hamiltonian which can only be 
written as a linear superposition of the Aa- Last expres- 
sion emphasizes that, while the states remain in the Sq 
space during the stochastic process, the average evolution 
can now simulate the evolution with an effective Hamil- 
tonian containing not only linear but also quadratic in 

Aa- 

The goal is now to take advantage of this generaliza- 
tion and reduce further the distance between the average 
evolution and the exact one. The most natural gener- 
alization of the last section is to minimize the average 
action 



1. Effective Hamiltonian dynamics deduced from the 
minimization 

The variations with respect to and Sdap give two 
sets of coupled equations between dq'^ and dtOap- The 
formal solution of the minimization can however be ob- 
tained by making an appropriate change on the varia- 
tional parameters prior to the minimization. In the fol- 
lowing, we introduce the notation = AaAp + ApAa 
where v denotes (a, /3) with a < (3. Starting from the 
general form of the effective evolution H23|l , we dissociate 
the part which contributes to the evolution of the {Aa) 
from the rest. This could be done by introducing the 
projection operator Vi- Equation H23|) then reads 

- |^dzS^„-i-^dt^.(i-7'i)s,| |g425) 

where the new set of parameters dz"^ is given by 

dzl ^dql+Y, di^.Mj {Q, \ApB,\ Qa) . (26) 

Similarly, the average evolution (dQbl transforms into 

= {Qb\\^dz''c,*Aa+Y,dcTMl-Vi)^V) 

where dz^ is given by 

dz^* = dq''^* + d(T. {Qb \B,Ap\ Qa) M-J. (28) 

In the following, we write B'^ — {1 — Vi)B^ and — 
B„{\ — Vi). The great interest of this transformation is 
to have {AaB^} = and {B'^Aa) = for all a and v. 
Accordingly, the variations with respect to Sz^ and Sz^ 
lead to 



ih{Qb\Aa\dQa) = {Qb\AaH\Qa) 

(29) 

lh{dQb\Aa\Qa) = {Qb\HAa\Qa), 

which gives closed equations for the variations dz^ and 
dz^ which are decoupled from the evolution of dui^, and 
dcTjy . These equations are identical to the ones derived in 
step 1 and can be again inverted as 



Y.dz:Aa\Qa) = §-VlH\Q,,), 



S = Tr {{ihdf - ihdf -H}D), 



(24) 



{QblY^dZaAa 



dt 



{Qb\HVi 



(30) 
(31) 
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On the other hand, the variations with respect to ScFi, 
and Suj,, lead to 



dt 



ih{Qb\B'l\dQa) = {Qb\B'lH\Qa) , 



(32) 



ih{dQb\B'^\Qa) = {Qb\HB'^\Qa) 



which again gives closed equations for duj^, and da^. 
These equations can be formally solved by introducing 
the two projectors V2 and associated respectively to 
the subspaces of operators By{l — V\) and {l — Vi)B^. V2 
differs from V2 due to the fact that Bi, operators and Aa 
operators do not a priori commute. Then, the effective 
evolution given by eq. H23|) becomes 



In summary, using the additional parameters associated 
with the stochastic contribution as variational parame- 
ters for the average action given by eq. (|24|l . one can 
further reduce the distance between the simulated evo- 
lution and the exact solution. When gaussian noises are 
used, this is equivalent to impose that the evolution of 
the correlations between operators Aa obtained by aver- 
aging over different stochastic trajectories also matches 
the exact evolution. 

C. Step 3: Generalization 



\dQa) 



while 




Y,dzlAa + {l 



V2)H\Qa) 



Vi)Y,dLU,B, ) \Qa) 



{dQb 



in 



V'2). 



(33) 



(34) 



In both cases, the first part corresponds to the projection 
of the exact dynamics on the space of observable {Aa). 
The second term corresponds to the projection on the 
subspace of the observable {AaAp) "orthogonal" to the 
space of the {Aa). 



2. Interpretation in terms of observable evolution 

The variation with respect to an enlarged set of pa- 
rameters does a priori completely determine the deter- 
ministic and stochastic evolution of the two trial state 
vectors. The associated average Schroedinger evolution 
corresponds to a projected dynamics. The interpretation 
of the solution obtained by variational principle is rather 
clear in terms of observable evolution. Indeed, from the 
two variational conditions, we can easily deduce that 



d{Aa) 



[Aa.H]), 



If the Hamiltonian H applied to the trial state can be 
written as a quadratic Hamiltonian in terms of Aa and 
if the trial states form an overcomplete basis of the total 
Hilbert space, then the above procedure can provide with 
an exact stochastic reformulation of the problem. If it is 
not the case, the above methods can be generalized by in- 
troducing higher order stochastic variables. Considering 
now the more general form 



5qS^ = ka + -^Ca ' + -^Ca 



,[3] 



2] 



J3] 



we suppose now that the only non vanishing moments 
for d^a^ and dr]a^ are the moments of order k (which are 
then assumed to be proportional to dt). For instance, we 

assume that ^ ' verifies 



(35) 
(36) 



Then without going into details, we can easily general- 
ize the method presented in step 2 and deduce that the 
average evolutions of the trial states will be given by 



\dQa) - -^{Vl+V2+V3 + ■■■}H\Qa) 

in 

m = -^{Qb\H{ri + r'^ + v', + ---} 

where the first terms contain all the information on the evolution of the {Aa), the second terms contain all the 
information on the evolution of the {AaAf^) which is not accounted for by the first term, the third terms contain all 
the information on the evolution of the {AaApAj) which is not contained in the first two terms, ... The procedure 
described here gives an exact Monte-Carlo formulation of a given problem if the Hamiltonian H applied on \Qa) or 
{Qb\ can be written as a polynomial of Aa. If the polynomial is of order k, then the sum stops at Vk. 



D. Summary and discussion on practical 
applications 

In this section, we propose a systematic method to re- 
place a general quantum problem by stochastic processes 



within a restricted class of trial state vectors associated 
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to a set of observable Ac,. Using variational techniques, 
we show that a hierarchy of stochastic approximations 
can be obtained. This method insures that at the level 
k of the hierarchy, all moments of order k or below of 
the observable Aa evolve according to the exact Ehren- 
fest equation over a single time step. The Monte-Carlo 
formulation might becomes exact if the Hamiltonian ap- 
plied to the trial state writes as a polynomial of the Aa 
operators. 

Aside of the use of variational techniques, we end up 
with the following important conclusion: Given an ini- 
tial density D = \Qa) {Qb\ where both states belongs to 
a given class of trial states associated to a set of opera- 
tors Aa, we can always find a Monte- Carlo process which 
preserves the specific form of D and insures that expec- 
tations values of all moments of the Aa up to a certain 
order k evolve in average according to the Ehrenfest the- 
orem associated to the exact Hamiltonian at each time 
step and along each trajectory. 

This statement will be referred to as the "existence 
theorem" in the following. Such a general statement 
is very useful in practice to obtain stochastic processes. 
Indeed, the use of variational techniques might become 
rather complicated due to the large number of degrees 
of freedom involved. An alternative method is to take 
advantage of the natural link made between the average 
effective evolution deduced from the stochastic evolution 
and the phase-space dynamics. Indeed, according to the 
existence theorem, we know that at a given level k of 
approximation, the dynamics of each trial state can be 
simulated by an average effective Hamiltonian insuring 
that all moments of order k or below matches the exact 
evolution. In general, it is easier to express the exact 
evolution of the moments and then 'guess' the associated 
SSE. 

The second important remark is that the above method 
can also be used to provide a SSE which maintains 
Tr{D) = 1 along each path. In that case, the action H24() 
can still be used but the density D should be replaced by 



D = 



\Qa) {Qb\ 
{Qb\Qa) 



(37) 



All the equations given above can then be equivalently 
derived. The main difference being that Aa is now re- 
placed by {Aa — (Aa)) in the variations of the states 
(equation (|14|l ). As we will see, the possibility to use 
constant trace formulation usually simplifies the use of 
the Ehrenfest theorem to guess the involved stochastic 
equations. 



phase-space evolution and the associated SSE is pointed 
out. 



A. The free wave-packet expansion 

We first consider a system initially in the ground state 
of a harmonic oscillator described by the Hamiltonian 
Ho ~ hu) (a+a -f |) where a and a+ are the usual cre- 
ation/annihilation operators with [a, a+] = 1. Denot- 
ing by 77 = muj/h, the initial wave-packet corresponds 
to a gaussian wave- function of width < >— At 
time i = 0, the harmonic constraint is released. During 
the free exp ansion, the wave-packet remains gaussian and 
verifies Hsl 



where m is the mass of the system. 



1. Stochastic formulation 

We show that this free expansion can be simulated 
by quantum diffusion between gaussian wave-functions of 
fixed widths. Let us assume that, at a given intermediate 
time step of the stochastic process, the density is given 
by averaging over densities of the form (which includes 
the initial condition) 



D = 



(/3|a)' 



(39) 



where \a) and (/?| are Bargmann states of the initial os- 
cillator defined as |q;) = e"°^ |0)[1|- These states verify 
a I a) = a la) and (/?| a+ = (/3|/3*. We use Bargmann 
states instead of coherent states to simplify the discus- 
sion. In particular, we have the transformation proper- 
ties 

|a-f da) = e'^""'' |a) , (/3 -I- = (/?| e'*^*", (40) 

which illustrates that the generator of transformations 
between Bargmann states are the a and a+ operators. 

Once the harmonic constraint is released the Hamilto- 
nian becomes the free Hamiltonian 



(41) 



IV. APPLICATIONS 



The action of H on \a) and {(3\ can be written respec- 
tively as 



In the following, we give few examples which have been 
studied recently in the context of Monte-Carlo techniques 
to illustrate how formal results obtained in previous sec- 
tion can be applied. In each case, the link between the 



H\a) 



-^(a+^ - 2aa^ + a^) \a) , 



H - -^(/3|(a2-2/3*a-t-/?*^), 
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which corresponds to polynomial of order 2 in a"*" and 
a respectively. Therefore, we expect that the exact dy- 
namics can be simulated with gaussian stochastic pro- 
cess between densities given by eq. (|39|) . The associated 
stochastic evolution of a and (3* are written as 

( da = d^ + d£,^^\ 
{ dp* = dl3*+d4'^\ 

where da and dj3* correspond to the deterministic part 
of the evolution, while c?^'^^ and dyyl^l are two indepen- 
dent gaussian stochastic variables with zero mean and 
following Ito rules 42]. In order to precise the nature 
of the stochastic evolution, one can use the fact that 
the first and second moments of the a and a+ opera- 
tors follow in average the exact evolution over one time 
step. Noting that, for D given by eq. H39() . we have 
(a+) = Tr{a+D) = /3* and (a) = Tr{aD) = a, the 
Ehrenfest theorem applied to (a+) and (a) gives the two 
equations 

d^ = dp^ =i^m* -a). (42) 
2m 

The application of the Ehrenfest theorem respectively to 
( ), ( ) and (a+a) gives 



d(a2) = 


m \ 




— a)a -1 




d{(3*^) = 


'^dtl 
m \ 




- a)P* 


-1} 


d{aP*) = 


'^dt{{P*^ 
m 


- a){a - 





Last condition is automatically fulfilled if we assume that 
d^I^l and d?]^^^ are independent stochastic variables. Ac- 
cording to Ito rules, we also have 

d{^ = 2ad^ + d(^^^d^^'^\ 

= 2l3*dj¥ + dri^^Ur]^^\ 

from which we deduce 

^^[2] ^^[2] ^ _rf^[2]rf^[2] ^ (43) 

2m 

A stochastic evolution of a and /3* compatible with the 
above conditions is then 




where dWi and dW2 are independent stochastic variables 
with zero mean and dWidWi = -dW2dW2 = f . The 
above c- number dynamic defines a stochastic evolution of 
the trial states through the relation H4U|) or equivalently 
can be interpreted as a diffusion process in the space of 
densities given by expression (|39|l . 



2. Recovering the exact dynamics from the average 
evolution 

Let us now check that the stochastic dynamics properly 
describes the exact evolution. For this, it is convenient 
to consider the stochastic complex variables defined as 
X = Tr{Dx) and P — Tr{Dp). For densities given by 
eq. ll^ . X and P are related to a and /3* through the 
relations 

P = ihyqiP* -a). 

Accordingly, X and P are stochastic complex variables. 
Starting from eq. (|45|l . we obtain 

dX = ^dt + dxu 

(46) 

dP = dX2, 

where the new stochastic variables dxi and dx2 verify 

dxi = dx2 ^ 0, 
dxidxi = dx2dx2 = 0, 

dXidX2 = T^—dt. 

Zm 

The Langevin equations (|46|l identifies with the classical 
equation of motion of the free particle with extra complex 
gaussian noises. Assuming that initially X{0) — P{0) — 
0, and using Ito rules of stochastic calculus, we deduce 
that 

Since along each stochastic path Tr{Dx^) = ^ +X'^, we 
finally end up with 

TKD^) = ^ + f^2t', 
Zi] Zm^ 

which is nothing but the exact result given by eq. (|38|l . 

With this simple example, we have illustrated how 
a quantum problem can be transformed into a Monte- 
Carlo process between densities formed of two Bargmann 
states. A remarkable aspect of the method used here is 
that at no time it is needed to minimize directly the ac- 
tion. We only took advantage of both the existence the- 
orem and Ehrenfest evolution to obtain the Stochastic 
Schroedinger evolution of trial states. It should however 
be kept in mind that the above stochastic process is inti- 
mately connected to the projection technique described 
in previous section. 

B. The Kerr oscillator 

As a second example, we consider the model Hamilto- 
nian known as the Kerr oscillator 

H = huia+a + ^hK{a+fa^. (47) 
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Using either the positive-P representation |8| or more re- 
cent works based on quantum-field theoretical techniques 
01 , it can be shown that the quantum problem of a sys- 
tem evolving with H can be replaced by an ensemble of 
stochastic c-number evolutions. Let us follow the same 
strategy as in previous example, and assume that at a 
given time step the density can be recovered by averaging 
over densities given by eq. H39|l . Using similar arguments 
as before, we know that the exact Hamiltonian dynamics 
can be simulated by a gaussian diffusion process on a and 
/?*, which have the general form given by eq. H42|) . The 
Ehrenfest theorem applied on (a) and (a+) gives directly 

'dbt = -idt {lo + K [P* a]) a, 

while the evolutions of ( a'^ ), ( ) and (a+a) lead re- 
spectively to 

A? = -2idt{uj + K[(3*a]}a'^ -iKdta^, 
df3*^ = 2idt{Lu + K[f3*a]}P*^ + iKdt(3'^, 
= 0. 



d{aP 



From the two sets of equations, we deduce 



2] 



iKdta^, 



straightforward compared to ref. [44|. Note that above 
stochastic equations are non-linear stochastic differential 
equation which might be numerically difficult to handle 



equa' 



C. The generalized Kerr Hamiltonian 

Up to now, we have presented examples where the ex- 
act Hamiltonian can be recovered using gaussian noises. 
To illustrate further the simplicity of the method, we 
consider the following Hamiltonian 



H = hiua 



1 



1 



d?7[2]dj^[2] = iKdtf3*\ 



-hK^{a+fa' + -hK2{a+fa'', (48) 

which will be referred hereafter as the generalized Kerr 
Hamiltonian. According to the existence theorem, in or- 
der to simulate exactly the Hamiltonian, it is necessary to 
include an extra contribution to the c-number evolution 

da = rf^+rfe'^' +di^^\ 
d(3* = dp*+dr]^^^ +dr)^^\ 

where first and second moments of d^^^l and drj^"^^ cancel 
out while third moments are proportional to dt. As in 
previous section, the use of the Ehrenfest theorem for the 
first and second moments leads respectively to 



Therefore, a stochastic Monte-Carlo process compatible 
with the above condition is given by 

da = -idt {uj + aK [(3* a]} + VKha dWi 
dp* = idt {lu + K [P*a]} /3* + P*VKhdW2 

where again dWi and dW2 are independent stochastic 
variables with zero mean and dWidWi — —dW2dW2 — 
j^. The above stochastic process is exactly the same 
as the one derived in ref. [4J| using a completely differ- 
ent technique. Indeed, although the presented framework 
uses phase-space arguments, it has a strong connection 
with the Hamiltonian dynamics. A great advantage here 
is that the guess of the Monte-Carlo process is rather 



da = -idtiuj + Ki[P*a] + ^K2[P*af ^■a, 
W ^ idt\uj + Ki [P*a] + ^K2 [P*af } P* , 



and 



dal^ldal^l = -idta^iKi+ K2[P*a]), 
dp*^^Up*^^^ = idtp*^{Ki+K2[P*a]). 



In order to determine the extra contribution, one should 
use the third moments. The evolution of ( a^^ ) gives for 
instance 



J 



d{P*f = ip*U3 



to + A'l [P*a] + \k2 [P*af 



?,{Ki+K2 [P*a]) + K2 



Using the identity 



d{P*Y = 3P*^dp* + 3P*dP*^^^dP*^^^ + dp*^^^dp*^^^dP*^^\ 

I 



gives 

d/3*'^ld/3*f^ld/3*'^l = idtfiK2P*^. 



Similarly, the evolution of ( a'^ ) gives 

dal^^dal^ldal^l = -idthK2a^. 
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Therefore, the exact evolution of the generaUzed Kerr the stochastic process on the c-numbers 
oscillator can be simulated with coherent states using 



da = ~idt {uj + Ki [/3*a] + ^i^a [P*af )■ a + {d^2 + ^6) a, 
dp* =idt\Lu + Ki [ffa] + ii^2 [13* af \ (3* + (d?72 + d^ys) P* 



r 



where the d^; and dr^i are independent stochastic vari- 
ables verifying 

d£,i = drji = 0, 



d^sd^s = dr]3dri3 = 0, 

d^^ = -~d^^2^-idt{Ki + K2[(3*a]), 
d^ad^^d^s = -drisdrisdrji = -idtK2. 

The exactness of the above stochastic process can be 
checked using the fact that for any density given by equa- 
tion H39|l and the stochastic evolution of a and (3* given 
above, we have the average relation 

gda(a+-(a+»£, + J^gd/?* (a- (a) ) = ^ , D] , (49) 

where H is the total Hamiltonian given by eq. H48|l . This 
relation shows that the average evolution identifies with 
the exact evolution for a given D and a single time step. 
Note that here the use of (a+ — (a+)) and (a — (a)) is 



related to the constant trace, 
each paths. 



Tr{D) = 1, imposed along 



such a stochastic reformulation. We show here that the 
proposed method also leads to an exact reformulation of 
the exact problem in terms of jumps between densities 



D = \N: ^a) {N : ^ 



(51) 



where both \N : (pa) and : 0b) stand for two Bose- 
Einstein condensates wave-function and where Tr[D) = 
1 along each path. 



1. Phase-space stochastic dynamics 

Let us assume that the system is at initial time to in a- 
specific Bose condensate wave- function |iV : 0), where \4>) 
is the associated single-particle state. Due to the nature 
of the initial state, all the information is contained in the 
one-body degrees of freedom. It is then convenient to 
define the one-body density matrix {i\pi\i) = (^afaj'^, 
whose matrix elements can be written as 



Pi 



N\ 



(52) 



D. Interacting Boson systems 

In recent years, several works have demonstrated that 
bosons or fermions interacting through a two- 
body potential can be simulated exactly by a Monte- 
Carlo method using simple trial state vector. Let us 
first consider the general case of N interacting bosons 
described by the general two-body Hamiltonian 



\T\j)a^ 



a, + - 



\ X! 1^12 1 Ik) afa+aiakl50) 

ijkl 



where the and af correspond to a complete set of 
creation/annihilation operators which obey bosons com- 
mutation rules. It has been shown in ref. p^ . that 
Bose-Einstein condensates wave-functions can be used 
as trial wave-functions. In that case, the generator of 
the transformations between two non-orthogonal wave- 
function are the single-particle operators afaj. Since 
H can be written as a quadratic polynomial in terms 
of these operators, we can have guessed the existence of 



being the projector on the single particle state {(p). 
According to the existence theorem, wc know that an ex- 
act Monte-Carlo reformulation can be obtained using the 
first and second moments of a complete set of one-body 
operators. This is equivalent to express the evolution of 
the one and two-body densities over one time step. In- 
deed, the two-body density denoted by pi2 is defined as 
usual as {ij\pi2\kl) — {alaf ajUi) . For a pure Bose- 
Einstein condensate, pi2 reads 



Pi2 = iV(iV-l)P4l)P42), 



(53) 



where Ptj,{i) means that the projection is made on the 
i*^ particle. The evolution of pi and pi2 over one 
time step are given by the first two equations of the 
Bogolyubov-Born-Green-Kirwood-Yvon (BBGKY) hier- 
archy I47I li^ , which reduces in the present case to 
the two equations of motion 



= [ft.MF(l),Pl] 



(54) 



«^^Pl2 = [ft-MF(l) + ft-MF (2), P12] +Pl2- (55) 
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hMF denotes the mean-field Hamiltonian given by 

hMF =T+iN -l)Tr2{vi2P4'2)), (56) 

where Tr2 corresponds to the partial trace on the sec- 
ond particle. B12 is the extra contribution which treats 
correlations beyond mean-field and reads 

B12 = (1-P0(1))(1-P42))«12P12 

-Pi2^'i2(l-P0(l))(l-P0(2)). (57) 

Due to this extra contribution, the exact many-body den- 
sity differs from a pure Bose-Einstein wave-function after 
one time step evolution. 

Starting from the initial densities given by eq. H52() 
and (|53|l . it is possible to guess a phase-space stochastic 



dynamics that simulates exactly the evolution of pi and 
P12. This is indeed the case if we assume that evolves 
according to a stochastic equation that verifies 



dP^ = — [hMF,P4,] , 

dP^il)dP^{2) = ^(l-P4l))(l-P42))t.i2P41)P42) 
-P4l)P^{2)v,2{l - P41))(1 - n(2))- 

Following ref. [IM we can write V12 as a sum 
of squares of hermitian one-body operators: V12 — 
Oa(1)Oa(2). We see that a proper choice for the 
stochastic evolution of P^ is 



dP<s> ^ ^ [hMF,P,p 
in 



^rf^i'l (I -P^)OaP^ 



5]<1p^Oa(1-P0), 



(58) 



where dS}-^ and cZr^^ are two sets of independent stochas- 
tic variables with mean zero and verifying d^^^'^d^^} = 
^AV§and<l<l = -<5Avf 

2. Equivalent stochastic process in Hilbert space 

The above phase-space stochastic dynamics can be 
simulated through a stochastic Schroedinger equation 
on single-particle wave-functions. Then, P(p{t) = |(/)) (0| 
transforms into a set of couples of wave-packets P^ {t + 
dt) = + d(f)a) {(j) + d(j)i,\ = l^a) (^bl, where states 
evolves over one time step according to 

d\<i,a) = (§W+Ea4"(1- 10a) 

d{M = {(t>b\{~%hMF + T.xdrll^Ox{l-P^)). 

where |(/'a(io)) = |06(^o)) = 10)- It turns out that the 
above stochastic process preserves P| = P^ for each tra- 
jectory, i.e. {(j)a I (^b) = 1. 

This stochastic evolution corresponds to quantum 
jumps in many-body space which transform the initial 
pure Bose-Einstein condensate into an ensemble of den- 
sities given by 

|iV:0,)(7V:4|, (59) 

with (0a I (pb) — 1 (which implies Tr{D) ~ 1) and is a pri- 
ori valid only for the first time step. Let us now assume 
that we start from one of this many-body densities. It 
could be shown without difficulty that all the above equa- 
tion remains valid except that now P^ = \4>a) {4>b\- There- 
fore, the single-particle SSE given above can be used to 



describe the long time evolution of the system. Note fi- 
nally that the Monte-Carlo equations presented here dif- 
fers from the "constant trace" scheme described in ref. 
|22| in particular due to the use of a common mean-field 
Hamiltonian in the evolution of |0a) and 

E. Interacting Fermions systems 

We consider now a system of fermions interacting 
through a two-body Hamiltonian 

H = ^PyO+aj + ^ X! l^i^l kl) ajd^aiak, (60) 

ij ijkl 

where the creation/annihilation operators now obey 
fermionic commutation rules and where V12 is antisym- 
metric. Since the derivation of a Monte-Carlo theory is 
very similar to the boson case, we only give here the im- 
portant steps. 

Let us assume that at a given time step, the exact 
density can be recovered by averaging over an ensemble 
of densities 

D=\^a){^b\, (61) 

where both states correspond to Hartree-Fock states. If 
we denote by {|A)}i=i,Ar a-^d {|ai)}^^^ the set of N 
single-particle states, we assume in addition that for each 
couples of Hartree-Fock states, associated singles-particle 
wave-functions verify {Pj \ai) = Sij. Accordingly, the 
one-body density matrix associated to a given D reads 

pi=5]|a,)(A|. (62) 
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We can easily verify that pf = pi and Tr{D) = 1. For 
each D given by eq. (|61|l . the two-body density is given 
by Pi2 = (1 — Pi2)piP2 where P12 corresponds to the 



permutation operator. The evolution of pi and pi2 over 
one time step are given by the two first equations of the 
BBGKY hierarchy which reads in that case 



— [hMF,Pl], 

= [hMF{l) + /imf(2), P12] 

+ (1 - pi)(l - p2)vi2PlP2 - PlP2Vl2{l - Pl){l - P2), 
I 



(63) 
(64) 



where the mean- field Hamiltonian now reads Hmf — 
T + Tr2{vi2P2)- This expression identifies with the one 
obtained for bosons except that is now replaced by the 
one-body density. Again, introducing a complete set of 
hermitian operators Oa, the previous expression can be 
simulated by a stochastic dynamics in phase-space given 
by eq. (where again is replaced by pi) equivalent 
to a stochastic process for single-particle wave-functions 
given by 

dm = mi^-^hMF+Y^dri^^ox{i~pi)^. 

This stochastic equation on single-particle wave- 
functions, which has been found using a different tech- 
nique in ref. 5lj, preserves the property {f3j \ ai) — 
dij. Therefore, it corresponds in many-body space to a 
Monte-Carlo procedure which transforms the initial set 
of densities into another set of densities with identical 
properties. 

1. Possible Extensions for Many-Body systems 

In the context of many-body systems, the above deriva- 
tions are restricted to particles interacting through two- 
body interactions. The method presented here can also 
be used to provide Monte-Carlo formulations when par- 
ticles interact through three-body or more interactions. 
Then higher order noises are necessary and evolution of 
the three-body or more density matrices should be esti- 
mated along stochastic paths. 

Another possible application of the method is to in- 
clude pairing correlation in the trial states. In that case, 
the relevant degrees of freedom along the path are not 
only the one-body density pij — Tr{afajD) but also the 
abnormal density defined as = Tr{ajaiD) (see for 
instance [solls^). Although we do not develop here ex- 
plicitly the formalism, it would be interesting to make 
the connection between the technique described above 



and the recent related works which use either directly 
a stochastic simulation in pha se-space or stochastic 
Schroedinger equations [53Ll54j . 

We gues s that the result is the one obtained recently 
in ref. [5J|. In that case, it has been shown that the ex- 
act two-body problem can be mapped into a quantum 
Monte-Carlo process between densities D — {^b\ 
written as a product of two Hartree-Fock Bogolyubov 
states with | ^a) = 1 along the path. 



V. DISCUSSION 

Different examples given above illustrate how one can 
take advantage of both phase-space evolution and Monte- 
Carlo process in Hilbert space. Most often, the ex- 
act problem is replaced by a set of coupled non-linear 
stochastic equations. Therefore, we would like to men- 
tion that we can end with same difficulties encountered 
sometimes in that case 8] with the appearance of unsta- 
ble trajectories. This is the case for instance with the 
stochastic nonlinear equations derived in section (|IVBp . 
which are numerically difficult to integrate j45l |. 

Another example is given by the two-mode Hamilto- 
nian of ref. j55j given by 



H 



(c+C2 + c+ci) + hK [ic+)^cl + {c+fc: 



where (c^, ci) and (cj, C2) correspond respectively to the 
creation/annihilation operators of two orthogonal sin- 
gle particles states denoted by \ui) and \u2), and obey 
boson commutation rules. This Hamiltonian has been 
used to mimic the dynamics of two coupled conden- 
sates. The exact static and dynamical solutions asso- 
ciated to this Hamiltonian can be obtained exactly and 
it has been retained as a test case in ref. [23| for Monte- 
Carlo methods. The framework presented here leads 
to an alternative stochastic processes between densities 
D = \N : (pa) {N : <j)i,\, where single-particle states \(j)a) 
and \(j)b) can be decomposed as 

|0a) = aa{t)\ui) + f3a{t)\u2) , 

\<l>b) = ab{t)\ui) + (3b{t)\u2) ■ 
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A direct application of the formalism developed for \/?ii^cJc2 leads to stochastic equations for and (3a 
bosons in section llV Ul with Oi — 



^hKcfci and O2 



daa 
dfia 



2 

idtQ 



Pa - 2idtK{N - 1) [alaa] Ua + ^/2Kh d^i ' - dW 



Ya - 2ldtK{N - 1) [PlPa] Pa + ^2^ 



dC2 



dW 



Oia, 
Pa. 



r 



where dW 



dCf^alaa+d^^2^P*,pa. 



These equations have 



to be completed with similar equations for al and P^ . It 
could in particular be checked that Tr{P^) — ct'lcta+PtPa 
is constant along each path. 

An illustration of the result obtain by a direct numer- 
ical integration of these equation is shown in figure ^ for 
the special case where all bosons are initially in the same 
state 1^2) and for the same parameters values taken as 
in ref. i.e. iV = 17, = 0.11^ and Qdt = 10^. 

Figure n presents the evolution of pi = (c^ci) as a func- 
tion of time. The stochastic calculation does perfectly 



chanics both with Stochastic Schroedinger Equation |23| 
or stochastic evolution in phase-space Therefore, to 
take fully advantage of these techniques one should de- 
velop specific numerical methods. This has been done 
for instance in refs. [2^ lisl Is^ using the fact that 
stochastic equations are generally not unique. Such a 
non-uniqueness also exists in the formalism described 
here. 



VI. CONCLUSION 




0.2 



0.15 



0.1 



0.05 



nt 

FIG. 1: Evolution of pi = ( c^ci ) as a function of time (cir- 
cles) obtained by averaging over 10^ stochastic trajectories. 
The result is compared to the exact evolution (solid lines). 
Error bars represent standard deviation. 

reproduce the exact results up to fit ~ 2.8 with a rather 
small number of trajectories. After this time, large de- 
viations with respect to the exact evolution occur, and 
could not be reduced by increasing the number of tra- 
jectories. A careful analysis show that some trajectories 
becomes hard to integrate numerically. Indeed, for some 
trajectories while Tr{P^) = 1, ||a^aa|| and ||/3J/?a|| may 
become very large. Accordingly, the implementation of 
the above stochastic equations requires the time step to 
be very small or, alternatively to develop specific numer- 
ical techniques. 

This last example illustrates that the stochastic equa- 
tions might be difficult to integrate numerically due to 
their non-linearity. This is a problem which seems to 
be recurrent in the context of quantum stochastic me- 



The main goal of this paper is to show that one can 
take advantage of the phase-space evolution to construct 
Monte-Carlo processes in Hilbert space of a restricted 
class of trial wave- functions. In the first part of this 
work, we show that given a class of trial wave-functions 
\Qa) associated to a set of specific operators Ao,, it is 
always possible to obtain a hierarchy of stochastic ap- 
proximations of a quantum problem in terms of quan- 
tum jumps between densities formed of couples of trial 
states D = \Qa) {Qb\- At the level k of the hierarchy, 
the existence of such a stochastic process is proved using 
variational techniques. The quantum diffusion obtained 
in such a way has a clear interpretation in phase-space 
evolution. Indeed, at a given level k we show that mo- 
ments of Aa calculated by averaging expectation values 
(^Q,j), {Aa^Aa^) {Aa^---Aa^) ovcr different trajec- 
tories match the exact evolution. The stochastic formu- 
lation can eventually be exact if the Hamiltonian applied 
to the trial state can be recast as a polynomial of Ac . 

The proof of the existence of such a hierarchy of 
stochastic formulation is very helpful to bridge stochas- 
tic mechanics in Hilbert space and phase-space evolution. 
In the second part of this article, several examples illus- 
trates how the Ehrenfest theorem can directly be used 
to guess stochastic Schroedinger equations. Finally, a 
critical discussion on numerical aspects is given. 
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